Super-Rough Glassy Phase of the Random Field XY Model in Two Dimensions 
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We study both analytically, using the renormalization group (RG) to two loop order, and numeri- 
cally, using an exact polynomial algorithm, the disorder-induced glass phase of the two-dimensional 
XY model with quenched random symmetry-breaking fields and without vortices. In the super- 
rough glassy phase, i.e. below the critical temperature T c , the disorder an d thermally averaged 
correlation function B(r) of the phase field #(x), B(r) = {[#(x) — #(x + r)] 2 ) behaves, for r 3> a, as 
B(r) ~ A{t) ln 2 (r/a) where r = |r| and a is a microscopic length scale. We derive the RG equations 
up to cubic order in r = (T c — T) /T c and predict the universal amplitude A(t) = 2r 2 — 2r 3 + C(r 4 ). 
The universality of A(t) results from nontrivial cancellations between nonuniversal constants of RG 
equations. Using an exact polynomial algorithm on an equivalent dimer version of the model we 
compute A(t) numerically and obtain a remarkable agreement with our analytical prediction, up to 
r^0.5. 



Disordered elastic systems are relevant to describe var- 
ious experimental situations ranging, for interfaces, from 
domain walls in ferromagnetic|l| or ferroelectric [2[ sys- 
tems, contact lines in wetting [Jto propagating cracks [H 
and, for periodic structures, from vortex lattices (VLs) 
in type-II superconductors [5| and Wigner crystals [6( , to 
charge or spin density waves Q. In most of these sys- 
tems, the large scale properties are described by a zero 
temperature fixed point, which can be described analyti- 
cally using the functional renormalization group . This 
latter has led to very accurate predictions, e.g. concern- 
ing various exponents, which could be, in some cases, 
successfully confronted to experiments or numerical sim- 
ulations [9(. 

In some cases, however, thermal fluctuations play an 
important role: this is the case for systems where the 
exponent describing the scale dependence of the free en- 
ergy fluctuations AF ~ L 6 is 8 = 0. It is then crucial to 
study the interplay between disorder and thermal fluc- 
tuations. While at zero temperature, Monte Carlo (MC) 
simulations, which are hampered by extremely long equi- 
libration times, could be circumvented by the use of pow- 
erful algorithms to compute directly the ground states 
using combinatorial optimization, the latter are of lit- 
tle use to study finite temperature properties. Here we 
consider a prototype of such situations, the classical 2D 
XY model with quenched random fields, known as the 
Cardy-Ostlund (CO) model (T^. It describes a wide class 
of systems including 2D periodic disordered elastic sys- 
tems, _such_as a randomly pinned planar array of vortex 
lines 
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_, LL5| , surfaces of crystals with quenched disorder 
[lij , random bond dimer models 03 and noninteracting 
disordered fermions in 2D [ll|, Il2j . In terms of a real 
phase field 0(x) G (—oo, oo), the CO model is defined by 
the partition function Z = J V9 er H l T with the Hamil- 
tonian 
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FIG. 1. The amplitude A(t), characterizing the super-rough 
phase. The squares indicate the numerical estimates obtained 
here using an exact polynomial algorithm [l(j. The 'one loop' 
curve indicates the one loop result A(r) = 2r 2 while the 'two 
loop' curve shows the two loop result © obtained here (we 
also show a Pade resummation of it). Aa is the result obtained 
in Ref. [Til from translating to (0 the free fermion calculation 
of Ref. [12j. We also show the values obtained at T = in 
the corresponding references. 



Here k is the elastic constant, a is the short-length-scale 
cutoff, and f and £ are quenched Gaussian random fields. 
Their nonzero correlations are given by 



/ l WP'(y) 



,2 9 



£(x)e(y)=T 2 ^-«5(x- y ), 



(2) 
(3) 



where i,j G {1,2} denote the components of f, T is 
the temperature, and . . . denotes the disorder aver- 
age. The disorder f must be introduced in the model 
as it is generated by the symmetry-breaking field under 



coarse graining [13J. The CO model exhibits a transi- 
tion at a critical temperature T c — A-kk between a high- 
temperature phase, where disorder is irrelevant, and a 
low-temperature disorder induced glass phase. It is de- 
scribed by a line of fixed points indexed by T, which, 
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thanks to the statistical tilt symmetry (STS), is not 
renormalized (9 = to all orders in perturbation theory). 
It displays many features of glassy systems, e.g., univer- 
sal susceptibility fluctuations 



IS] , and aging 



The most striking effect of disorder on the statics con- 
cerns the two-point correlation function (CF) B(r) — 
([0(x) - 6>(x + r)] 2 ). While for T > T c the interface is 
logarithmically rough, B(r) ~ AT/T C In (r/a), it becomes 
super rough for T < T c where 



B(r) = A(t) In 2 (r/a) + 0[\n(r/a)] , 



(4) 



with r = (T c — T)/T c . The temperature T is determined 
by the connected CF B c (r) 
(4T/T C ) In (r/a). 



;[0(x)-0(x + r)P) c ~ 
A physical realization of (JTJ is the 
VL confined in a superconducting film with a parallel 
magnetic field [3, [5l|, [22j]. Such a geometry was re- 



alized experimentally on mesoscopic devices 15]. The 
ln 2 (r) growth in (0} results in a loss of translational 
order. Indeed it can be shown [l(| [23[ that d3| im- 
plies that the CF of the order parameter for transla- 
tional order of the VL decays faster than any power law 
In ^ e ig(0(x+r)-e(x))^ ^ ln 2 (r). This can be probed by 

neutron scattering experiments or direct observation via 
STM or scanning superconducting quantum interference 
device probes [241 ] . 

Although the amplitude A(t) of this intri guing ln 2 (r ) 
has been the subject of numerous studies [18t 
none of them was able to establish a quantitative com- 
parison between analytical results and numerical simu- 
lations, and there are several reasons for this gap. The 
first one is that A(t) was, analytically, only known, at 
lowest order in r, A(t) — 2r 2 + 0(t 3 ) 29]: its domain of 
validity is thus restricted to a narrow region close to T c , 
where the amplitude of the ln 2 (r) term is small and thus 
hard to isolate accurately from the subleading logarith- 
mic correction in (|4]) . The second reason is that numerics 
is very delicate, given that standard MC simulations are 
quite inefficient for T < T c . Fortunately, there exists an 
exact polynomial algorithm, called the domino shuffling 
algorithm (DSA), which allows us to sample directly the 
related random bond dimer model, without running MC 
simulations. This algorithm was used in Ref. jl8j |. which 
showed that A(t) oc t 2 , without estimating the prefac- 
tor. Notice that, in its original formulation as used in 
Ref. [18|] , the DSA suffers from strong finite size effects 
reminiscent of the arctic circle phenomenon |30j in the 
pure dimer model. 

In this Letter, we perform a quantitative comparison, 
in a wide temperature range, between analytical and nu- 
merical predictions. Such a comparison is rendered pos- 
sible (i) thanks to a precise calculation of A(t), using 
various RG schemes, yielding the following expression to 
two loop order: 



A(t) = 2r 2 - 2r 3 + 0(t a 



(5) 



and (ii) thanks to a careful Fourier analysis of the two- 
point CF. It is computed here using an improvement of 
the DSA [3l|, where finite size effects are significantly 
reduced. The result of this comparison is shown in 
Fig. [TJ We see a remarkable agreement between both 
approaches, even far beyond T c down to r ~ 0.5. 

Our analytical study determines the RG equations for 
the model (TT]) to two-loop order. We use the replica 
method to treat the disorder [32| and obtain the repli- 
cated Hamiltonian H rep = H^ ep + if[ ep , with the har- 
monic part 
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E / & x {£f S <* [(W Q (x)) 2 +m 2 (e Q (x)) 2 ] 
_£v0 Q (x).V^(x)}. (6) 



The mass m is introduced as an infrared cutoff while a 
and j5 denote n replica indices. Both to and n are set to 
zero at the end. The anharmonic part reads 



rep 
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E J d 2 xcos(9 a (x) - 9 p (x)). (7) 

a/3 



We compute the two CFs: 



0(x) = (0(x)0(Q))h - (9(x)) H (9(0)) H , 



So(x) = <0(x)) ff <0(O)) ff , 



(8) 
(9) 



where the former measures the (disorder averaged) ther- 
mal fluctuations while the latter measures the fluctua- 
tions due to disorder of the (thermally averaged) phase 
field. These CFs can be obtained from CFs of repli- 
cated fields by decomposing G a /3(x) := ((6 a (x)9p(Q))) = 
8 a f3G{x) + Go( x ), where G(x) is called the connected part 
and C?o( x ) the off-diagonal part. To compute them we 
use the harmonic part Hq &p as the "free" theory and 
treat 7J[ ep in perturbation theory in g. Here we denote 
by ((..}) averages over the complete Hamiltonian H rep 
and by (..) averages over the free part H^ ep . 

We start by computing the CF for g = and find 



G a/ 3(q) 



T S, 



a/3 



aT 2 



K q 2 + to 2 2-7TK 2 (q 2 + TO 2 )' 



+ 0(n). (10) 



In real space one obtains G Q| g(x) = (0 Q (x)6' ( g(O)) = 
6 a pG(x) + Go(x). The connected part behaves at small 
distances |x| <C (cto) -1 as 



G(x) = -(1 - r) In [c 2 m 2 (a; 2 + a 2 )] 



(11) 



with c = e 7B /2 and is the Euler constant. In 
([TT]) we have introduced the ultraviolet regularization 
by the parameter a (33|. The off-diagonal part of 
the CF at small distances reads Gq(x) = — 2<r(l — 
t) 2 In [ec 2 TO 2 (x 2 + a 2 )] + 0(n). 

The STS of the model manifests itself by the in- 
variance of the non-linear part 7J[ e?) under the change 
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x) — »• 9 a (x) + (/>(x) for an arbitrary </>(x) . As discussed 
1J,|25|,|34| it implies two properties: (i) Gq(x) does not 



appear to any order in perturbation theory in g and (ii) 
the disorder-averaged thermal CF is <?(x) = G(x) to all 
orders in g. This implies that T can be measured from 
the amplitude of the logarithm in Q{x) ~ 2(1 — r) In x 
at large x. Because of property (i) Go( x ) only receives 
additive corrections, e.g. corrections to a which, in the 
present model, change its above logarithmic behavior into 
a squared- logarithm behavior for Qq(x), obtained below. 

To obtain the scaling equations beyond lowest order we 
compute the effective action up to 0(g 3 ), which reads 

r =E/ [W) 2 +- 2 (^) 2 ] 

-V6» Q • V0« - ^-c 2 m 2 cos((9 Q - (12) 
2ir > 



a/3 
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47T 

in terms of renormalized couplings g R and a R . Their ex- 
plicit dependence on the bare parameters leads to the fol- 
lowing scaling equations in terms of the scale £ = — him: 



dg R 
d£ 

d<J R 
d£ 



-j^ = 2t 9r - Ag% - Brg\ + Cg 3 R , 



Dgl + Erg R -Fg% 



(13) 
(14) 



and dr/d£ — 0, which encodes the exact result Q{x) — 
G(x) from STS. Remarkably, although A, B, . . . , F are 
nonuniversal constants, they satisfy the universal ratios 

A 2 /D = 8, A 2 /C = 4, F + BD - AE/2 = , (15) 

which will ultimately enter into the expression of phys- 
ical quantities, including A{t) in (UJ). We have ob- 
tained these values through three different regulariza- 
tion schemes, for details see Ref. [IH . These equations 
generalize to two-loop order the one-loop equations ob- 
tained in [H [H, M, From ([T3]) we see that 
the model has a transition at r = 0, i.e. T = T c . 
For T > T c the renormalized coupling g R (£) flows to 
zero, while for T < T c it flows to a finite value g R = 
2t/A + (4C - 2AB)t 2 /A 3 + 0(t 3 ). The asymptotic so- 
lution of (0) is a R (£) ~ cr + [da(g* R )/d£]£: while a 
is nonuniversal and leads only to logarithmic growth of 
the CF, the second term yields the ln 2 (r) growth [lij ]. 
To estimate the off-diagonal CF at a given wave vec- 
tor q, one considers the limit m q and argues that 
q sets the scale £* = ln[l/ (aq)] at which one stops the 
flow. Replacing a by its effective value at that scale, 
i.e. a — > <t r (£) we get, from (Ti"0| the small q behavior 
& (q) ~ Ml -T) 2 dai ^ ] ln ^ a >1 , which leads to © and 
([5]) . The amplitude ([5]) is universal thanks to the remark- 
able combination of nonuniversal constants in da(g R ) / d£. 
Equation (j4]) can be obtained more rigorously by calcu- 
lating the two-point function [35j ] . 

We have performed simulations to estimate numeri- 
cally the amplitude A(t) and compare it with (f5]). For 



that purpose, we use the mapping between (fTl and a 
weighted dimer model defined on a 2D lattice [H [3, 
for which there exists a polynomial DSA [lfj. For tech- 
nical reasons, it is designed for a special lattice Aj, called 
the Aztec diamond of size L [see Fig. Ufa)]. To each 
bond between nearest neighbors (r,r') on Al we assign 
a quenched random variable e r y : here we consider inde- 
pendent Gaussian variables of zero mean and unit vari- 
ance. The dimer model consists of all complete dimer 
coverings of Al, where the weight W(C) of a dimer cov- 
ering C is given by 

W(C) = exp (-H d /T d ) , H d = ]T e r , r , (16) 

L ^ e ' (r,r')GC 

where Z^ie) is the partition function. Hence the limit 
T d — > corresponds to a "strong" disorder regime while 
the limit T d — >■ oo corresponds to dimer coverings with 
uniform weights. The DSA generates uncorrelated dimer 
configurations, directly sampled with the equilibrium 
weight (HHJ), without the need to run a slow MC algo- 
rithm. In addition, this is a polynomial algorithm (with 
a computational time ~ L 3 ). The dimer covering of the 
Aztec diamond is, however, known to suffer from strong 
finite size effects [3fl |. Here we minimize significantly 
these effects by using a recent improvement of the DSA 
which allows for the existence of bonds with zero weight 
3l| . We use it to study the random dimer model directly 
on a square lattice, which exhibits less pronounced finite- 
size effects [37j. This algorithm is very flexible and will 
be very useful to study other dimer systems in various 
2D geometries. 

To a given dimer covering C, we assign a discrete height 
field, defined on the center of the squares [see Fig . [Ha)l, 
i.e. on the dual lattice A® of Al, as follows 38]. The 
bonds of A® are oriented such that the unit cells of A® 
that enclose the blue sites of Al are circled counterclock- 
wise. Assign —3 to the difference of neighboring heights 
along the oriented bonds if a dimer is crossed and +1 oth- 
erwise. This yields single-valued heights up to an over- 
all constant, the heights on the boundaries of Al being 
then fixed as in Fig. Ufa). This defines a height field 
H r = Hij, with r = m x + ju y and the relative height 
h = H— -< H >~ where -< H y is a spatial average of 
H over Al- For uniform dimer coverings, corresponding 
to e r r / = or T d — > oo, one can show that the fluc- 
tuations of h in the continuum limit (and in the bulk) 
are described by a Gaussian free field [H, [3§| , i. e. by 
the Hamiltonian in (fT]) without disorder (f = 0, £ = 0) 
at r = 0. For inhomogeneous random bonds e r r r one 
expects instead that in the continuum limit, the fluctu- 
ations of h are described by the CO model (fT]) with the 
substitution 6 ->• h x 2tt/4 [ijj [3. This factor 2tt/4 is 
required because the energy associated to the height con- 
figurations (fT6|) is invariant under a global shift h — > h+4. 

The temperature Td of the dimer model does not coin- 
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FIG. 2. a) Dimer covering of an Aztec diamond of size 4, A&. 
The blue points allow us to define the height field, which are 
the integers on the dual lattice Af. b) Plot of q 2 C(q) as a 
function of In q. The squares correspond to our numerical data 
obtained for lattice size L = 384. The slope of the straight 
line indexed by 'one loop' and 'two loop' is given respectively 
by the one loop A(t) = 2r 2 and the two loop estimate 
while As corresponds to Ref. [l^]. Without a In 2 term in Q, 
one would expect a straight line with vanishing slope: this is 
ruled out by our data and demonstrates the existence of the 
super-rough phase. 



cide with the temperature of ([TJ. To compute r we use 
the STS (fTTj) and measure 



(hr){hr) ~2(l-r)lni, (17) 



which provides a precise estimate of r. We have checked 
that our numerical estimate is in good agreement with 
the analytic results for r and for other thermodynamical 
observables obtained in Ref. 



17] 



We want to compute numerically the amplitude of 
the ln 2 (r) term in (|T]). Extracting this amplitude pre- 
cisely from B(r) is however difficult, since the sublead- 
ing corrections are of order O(lnr). The calculation is 
more accurate in Fourier space [26|, [S(j, defining A q = 
L- 2 X) r /i r e iq r . The CF C(q) of these Fourier compo- 
nents is expected, from ((4]), to behave for small q as 



C(q) = (/i q )(/i- q > * -A(t) 



8 „ v In (1/?) 



0(q~ 2 ) , (18) 



where q — |q|. In Fig. [UL), we plot of q 2 C(q) as a 
function of lng for r « 0.33 (T^ = 0.25). These data 
have been obtained for a system size L = 384 by aver- 
aging over 10 5 realizations of the random bonds e r y's. 
They support the expected behavior in (fT8|) for small 
q, q < 1: they are indeed well described by a straight 
line, q 2 C(q) — — 8A(r)/-7rln(7 + b . Note that the down- 
wards bending for the smallest g's is a finite size effect. 
In Fig. [UJb) we also show four different straight lines 
corresponding to different couples [^4(t), bo]. The line in- 
dexed by 'Best fit' corresponds to the best fit of these 
data by a straight line: the value of A(t) obtained in 
this way corresponds to the squares on Fig. [TJ In the 
three other cases, the slope of this straight line is evalu- 
ated from the one- and two-loop ([5} results respectively, 



while the straight line indexed by 'Ag' corresponds to the 
slope computed in [ll[ from the result in Ref. [12], with 
As(t) = 2t 2 (1 — t) 2 . In all cases the constant bo is a 
fitting parameter. One clearly sees that the two loop re- 
sult is a significant improvement over the one loop result 
and describes very well our numerical data. Clearly, Ag 
underestimates our data. 

Let us now discuss the numerical results for A(r) in 
Fig. [TJ As compared to Ref. [18J, here we can discuss a 
much broader range of values of r which extends deep 
into the glass phase. First we observe that our two loop 
result is in very good agreement with our numerics up 
to r w 0.5. In contrast, the curve Ag(r) is signifi- 
cantly smaller than our numerical values and can be ruled 
out. For smaller temperature, r > 0.5 the discrepancy 
between ([5]) and the numerical value increases, as ex- 
pected. In Fig. [TJ we have also quoted the numerical val- 
ues which were obtained independently at zero tempera- 
ture by an exact ground state calculation for the solid-on- 
solid model on a disordered substrate in [26, 27, This 
model is also described, in the continuum limit, by the 
model ([TJ. In particular, our data match smoothly with 
the most recent numerical estimate obtained in Ref. [i(| , 
yielding A(t = 0) = 0.39. We also show an estimate 
based on a one-loop functional RG calculation at T = 
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The fact that the two loop formula (0 vanishes 
at T = is of course an unphysical feature, which can 
be cured by considering various guesses or Pade (i.e. ra- 
tional functions in r) approximations which have the 
same expansion as ([5]) to order r 3 . One such formula 
A(t) = 2t 2 (1 — r/2) 2 is plotted in Fig. [TJ being simple, 
it also has a reasonable structure to correct the result of 
Ref. 0. 

In conclusion, we have obtained an accurate descrip- 
tion of the glassy phase of the Cardy-Ostlund model 
down to temperatures T c /2 < T, with excellent agree- 
ment for the amplitude of the square logarithm between 
theory and numerics. Understanding the glass phase be- 
low T c /2 is an important challenge for the future. 
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